library(nlme)
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
##
## lmList
library(ggplot2)
library(sjPlot)
library(AER)
## Loading required package: car
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.12.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:survival':
##
## kidney
## The following object is masked from 'package:lme4':
##
## ngrps
## The following object is masked from 'package:stats':
##
## ar
# read dataset
df <- read.csv("alpha_diversity_dispersal.csv", row.names = 1, check.names = 0)
# only keep richness
df <- df[, -c(2:7)]
# remove samples from day 0
df <- subset(df, !grepl('___', rownames(df)))
# select soil from sea water treatment
df <- df[df$Water == "sea",]
df <- df[, -4]
df$Soil <- factor(df$Soil, levels = c("0", "70"))
df$Frequency <- factor(df$Frequency)
df$box <- factor(paste(df$Soil, df$Frequency, df$replicates, sep = "_"))
str(df)
## 'data.frame': 144 obs. of 6 variables:
## $ richness : int 828 1075 1101 1100 1108 1124 792 900 727 1074 ...
## $ Soil : Factor w/ 2 levels "0","70": 1 1 1 1 1 1 1 1 1 1 ...
## $ Frequency : Factor w/ 4 levels "f1","f2","f3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : int 12 12 12 16 16 16 2 2 2 20 ...
## $ replicates: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
## $ box : Factor w/ 24 levels "0_f1_A","0_f1_B",..: 1 2 3 1 2 3 1 2 3 1 ...
head(df)
## richness Soil Frequency Day replicates box
## 0_f1_sea_12_A 828 0 f1 12 A 0_f1_A
## 0_f1_sea_12_B 1075 0 f1 12 B 0_f1_B
## 0_f1_sea_12_C 1101 0 f1 12 C 0_f1_C
## 0_f1_sea_16_A 1100 0 f1 16 A 0_f1_A
## 0_f1_sea_16_B 1108 0 f1 16 B 0_f1_B
## 0_f1_sea_16_C 1124 0 f1 16 C 0_f1_C
# histogram
hist(df$richness, main = "histogram of richness")
(f <- ggplot(df, aes(x = Day, y = richness, shape = Frequency, color = Soil)) +
geom_point(size = 2, alpha=0.9) +
labs(x = "Days", y = "Richness")+
theme_light())
(f <- ggplot(df, aes(x = Day, y = richness, shape = Frequency, color = Soil)) +
geom_point(size = 2, alpha=0.9) +
facet_grid(Soil ~ Frequency) +
labs(x = "Days", y = "Richness")+
theme_light())
(f <- ggplot(df, aes(x = factor(Day), y = richness, shape = Frequency, color = Soil)) +
geom_boxplot() +
geom_point(size = 2, alpha=0.9, position = position_jitterdodge()) +
labs(x = "Days", y = "Richness")+
theme_light())
# full model
glm0 <- glm(richness~Soil*Day*Frequency, family="poisson", data = df)
summary(glm0)
##
## Call:
## glm(formula = richness ~ Soil * Day * Frequency, family = "poisson",
## data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -17.3265 -1.4315 0.2237 2.4846 8.3313
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.686904 0.015134 441.859 < 2e-16 ***
## Soil70 0.340859 0.019967 17.071 < 2e-16 ***
## Day 0.017359 0.001186 14.637 < 2e-16 ***
## Frequencyf2 0.163617 0.020703 7.903 2.72e-15 ***
## Frequencyf3 0.131680 0.020845 6.317 2.67e-10 ***
## Frequencyf4 0.068546 0.021003 3.264 0.001100 **
## Soil70:Day -0.007455 0.001579 -4.721 2.34e-06 ***
## Soil70:Frequencyf2 -0.124594 0.027741 -4.491 7.08e-06 ***
## Soil70:Frequencyf3 -0.104908 0.027846 -3.767 0.000165 ***
## Soil70:Frequencyf4 0.122734 0.027641 4.440 8.99e-06 ***
## Day:Frequencyf2 -0.004736 0.001633 -2.901 0.003722 **
## Day:Frequencyf3 -0.004244 0.001643 -2.583 0.009783 **
## Day:Frequencyf4 0.001474 0.001643 0.898 0.369436
## Soil70:Day:Frequencyf2 -0.004220 0.002215 -1.905 0.056825 .
## Soil70:Day:Frequencyf3 -0.002340 0.002219 -1.055 0.291584
## Soil70:Day:Frequencyf4 -0.018084 0.002206 -8.196 2.48e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4746.1 on 143 degrees of freedom
## Residual deviance: 2489.2 on 128 degrees of freedom
## AIC: 3796.2
##
## Number of Fisher Scoring iterations: 4
plot(glm0)
# check the mean and var of Richness
mean(df$richness); var(df$richness)
## [1] 1131.972
## [1] 36103.89
# check overdispersion
dispersiontest(glm0)
##
## Overdispersion test
##
## data: glm0
## z = 5.0738, p-value = 1.95e-07
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 16.18844
# update full model
glm1 <- glm(richness ~ Soil*Day*Frequency, family = "quasipoisson", data = df)
summary(glm1)
##
## Call:
## glm(formula = richness ~ Soil * Day * Frequency, family = "quasipoisson",
## data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -17.3265 -1.4315 0.2237 2.4846 8.3313
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.686904 0.064583 103.539 < 2e-16 ***
## Soil70 0.340859 0.085211 4.000 0.000106 ***
## Day 0.017359 0.005061 3.430 0.000813 ***
## Frequencyf2 0.163617 0.088350 1.852 0.066343 .
## Frequencyf3 0.131680 0.088959 1.480 0.141267
## Frequencyf4 0.068546 0.089633 0.765 0.445835
## Soil70:Day -0.007455 0.006738 -1.106 0.270660
## Soil70:Frequencyf2 -0.124594 0.118385 -1.052 0.294577
## Soil70:Frequencyf3 -0.104908 0.118833 -0.883 0.378987
## Soil70:Frequencyf4 0.122734 0.117961 1.040 0.300087
## Day:Frequencyf2 -0.004736 0.006967 -0.680 0.497894
## Day:Frequencyf3 -0.004244 0.007011 -0.605 0.546012
## Day:Frequencyf4 0.001474 0.007010 0.210 0.833756
## Soil70:Day:Frequencyf2 -0.004220 0.009454 -0.446 0.656126
## Soil70:Day:Frequencyf3 -0.002340 0.009471 -0.247 0.805201
## Soil70:Day:Frequencyf4 -0.018084 0.009416 -1.921 0.057010 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 18.21199)
##
## Null deviance: 4746.1 on 143 degrees of freedom
## Residual deviance: 2489.2 on 128 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
# update model by removing non-significant factors
glm2 <- glm(richness ~ Soil * Day, family = "quasipoisson", data = df)
summary(glm2)
##
## Call:
## glm(formula = richness ~ Soil * Day, family = "quasipoisson",
## data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -18.6252 -1.6266 0.4184 2.8346 7.7032
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.779524 0.031647 214.226 < 2e-16 ***
## Soil70 0.313968 0.042300 7.422 1.02e-11 ***
## Day 0.015424 0.002493 6.186 6.39e-09 ***
## Soil70:Day -0.013585 0.003392 -4.005 1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 18.98598)
##
## Null deviance: 4746.1 on 143 degrees of freedom
## Residual deviance: 2851.6 on 140 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
# full model
lm0 <- lm(log(richness) ~ Soil * Day * Frequency, data=df)
plot(lm0)
plot(lm0$residuals~df$Day)
plot(lm0$residuals~df$Soil)
plot(lm0$residuals~df$Frequency)
summary(lm0)
##
## Call:
## lm(formula = log(richness) ~ Soil * Day * Frequency, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.58917 -0.03811 0.01901 0.07819 0.23875
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.6708410 0.0638063 104.548 < 2e-16 ***
## Soil70 0.3309207 0.0902357 3.667 0.000358 ***
## Day 0.0180250 0.0052567 3.429 0.000816 ***
## Frequencyf2 0.1808577 0.0902357 2.004 0.047150 *
## Frequencyf3 0.1341841 0.0902357 1.487 0.139464
## Frequencyf4 0.0813477 0.0902357 0.902 0.369015
## Soil70:Day -0.0069091 0.0074341 -0.929 0.354441
## Soil70:Frequencyf2 -0.1202914 0.1276126 -0.943 0.347646
## Soil70:Frequencyf3 -0.1041991 0.1276126 -0.817 0.415716
## Soil70:Frequencyf4 0.1322008 0.1276126 1.036 0.302177
## Day:Frequencyf2 -0.0056138 0.0074341 -0.755 0.451554
## Day:Frequencyf3 -0.0041081 0.0074341 -0.553 0.581495
## Day:Frequencyf4 0.0006347 0.0074341 0.085 0.932092
## Soil70:Day:Frequencyf2 -0.0057465 0.0105134 -0.547 0.585616
## Soil70:Day:Frequencyf3 -0.0025173 0.0105134 -0.239 0.811148
## Soil70:Day:Frequencyf4 -0.0185600 0.0105134 -1.765 0.079887 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.142 on 128 degrees of freedom
## Multiple R-squared: 0.4508, Adjusted R-squared: 0.3864
## F-statistic: 7.004 on 15 and 128 DF, p-value: 5.976e-11
# update full model
lm0 <- lm(log(richness) ~ Soil * Day, data=df)
plot(lm0)
summary(lm0)
##
## Call:
## lm(formula = log(richness) ~ Soil * Day, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64275 -0.03683 0.01994 0.09358 0.22086
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.769938 0.032648 207.360 < 2e-16 ***
## Soil70 0.307848 0.046171 6.668 5.55e-10 ***
## Day 0.015753 0.002690 5.857 3.22e-08 ***
## Soil70:Day -0.013615 0.003804 -3.579 0.000474 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1453 on 140 degrees of freedom
## Multiple R-squared: 0.3709, Adjusted R-squared: 0.3575
## F-statistic: 27.52 on 3 and 140 DF, p-value: 4.742e-14
# first approach ---
# full model
lmm0_1 <- lme(log(richness) ~ Soil * Day *Frequency, random=~1|box, data=df, correlation=corAR1(form = ~Day|box))
summary(lmm0_1)
## Linear mixed-effects model fit by REML
## Data: df
## AIC BIC logLik
## -22.5273 31.66128 30.26365
##
## Random effects:
## Formula: ~1 | box
## (Intercept) Residual
## StdDev: 3.786054e-06 0.1420281
##
## Correlation Structure: ARMA(1,0)
## Formula: ~Day | box
## Parameter estimate(s):
## Phi1
## 0
## Fixed effects: log(richness) ~ Soil * Day * Frequency
## Value Std.Error DF t-value p-value
## (Intercept) 6.670841 0.06380629 112 104.54833 0.0000
## Soil70 0.330921 0.09023572 16 3.66729 0.0021
## Day 0.018025 0.00525670 112 3.42895 0.0008
## Frequencyf2 0.180858 0.09023572 16 2.00428 0.0623
## Frequencyf3 0.134184 0.09023572 16 1.48704 0.1564
## Frequencyf4 0.081348 0.09023572 16 0.90150 0.3807
## Soil70:Day -0.006909 0.00743409 112 -0.92938 0.3547
## Soil70:Frequencyf2 -0.120291 0.12761258 16 -0.94263 0.3599
## Soil70:Frequencyf3 -0.104199 0.12761258 16 -0.81653 0.4262
## Soil70:Frequencyf4 0.132201 0.12761258 16 1.03595 0.3156
## Day:Frequencyf2 -0.005614 0.00743409 112 -0.75514 0.4518
## Day:Frequencyf3 -0.004108 0.00743409 112 -0.55261 0.5816
## Day:Frequencyf4 0.000635 0.00743409 112 0.08538 0.9321
## Soil70:Day:Frequencyf2 -0.005746 0.01051339 112 -0.54658 0.5858
## Soil70:Day:Frequencyf3 -0.002517 0.01051339 112 -0.23944 0.8112
## Soil70:Day:Frequencyf4 -0.018560 0.01051339 112 -1.76536 0.0802
## Correlation:
## (Intr) Soil70 Day Frqnc2 Frqnc3 Frqnc4 Sl70:D S70:F2
## Soil70 -0.707
## Day -0.851 0.602
## Frequencyf2 -0.707 0.500 0.602
## Frequencyf3 -0.707 0.500 0.602 0.500
## Frequencyf4 -0.707 0.500 0.602 0.500 0.500
## Soil70:Day 0.602 -0.851 -0.707 -0.426 -0.426 -0.426
## Soil70:Frequencyf2 0.500 -0.707 -0.426 -0.707 -0.354 -0.354 0.602
## Soil70:Frequencyf3 0.500 -0.707 -0.426 -0.354 -0.707 -0.354 0.602 0.500
## Soil70:Frequencyf4 0.500 -0.707 -0.426 -0.354 -0.354 -0.707 0.602 0.500
## Day:Frequencyf2 0.602 -0.426 -0.707 -0.851 -0.426 -0.426 0.500 0.602
## Day:Frequencyf3 0.602 -0.426 -0.707 -0.426 -0.851 -0.426 0.500 0.301
## Day:Frequencyf4 0.602 -0.426 -0.707 -0.426 -0.426 -0.851 0.500 0.301
## Soil70:Day:Frequencyf2 -0.426 0.602 0.500 0.602 0.301 0.301 -0.707 -0.851
## Soil70:Day:Frequencyf3 -0.426 0.602 0.500 0.301 0.602 0.301 -0.707 -0.426
## Soil70:Day:Frequencyf4 -0.426 0.602 0.500 0.301 0.301 0.602 -0.707 -0.426
## S70:F3 S70:F4 Dy:Fr2 Dy:Fr3 Dy:Fr4 S70:D:F2 S70:D:F3
## Soil70
## Day
## Frequencyf2
## Frequencyf3
## Frequencyf4
## Soil70:Day
## Soil70:Frequencyf2
## Soil70:Frequencyf3
## Soil70:Frequencyf4 0.500
## Day:Frequencyf2 0.301 0.301
## Day:Frequencyf3 0.602 0.301 0.500
## Day:Frequencyf4 0.301 0.602 0.500 0.500
## Soil70:Day:Frequencyf2 -0.426 -0.426 -0.707 -0.354 -0.354
## Soil70:Day:Frequencyf3 -0.851 -0.426 -0.354 -0.707 -0.354 0.500
## Soil70:Day:Frequencyf4 -0.426 -0.851 -0.354 -0.354 -0.707 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.1482469 -0.2683141 0.1338492 0.5505475 1.6809933
##
## Number of Observations: 144
## Number of Groups: 24
# update model
lmm1_1 <- lme(log(richness) ~ Soil * Day, random=~1|box, data=df, correlation=corAR1(form = ~Day|box))
summary(lmm1_1)
## Linear mixed-effects model fit by REML
## Data: df
## AIC BIC logLik
## -104.2061 -83.61463 59.10306
##
## Random effects:
## Formula: ~1 | box
## (Intercept) Residual
## StdDev: 1.176309e-05 0.1453448
##
## Correlation Structure: ARMA(1,0)
## Formula: ~Day | box
## Parameter estimate(s):
## Phi1
## 0
## Fixed effects: log(richness) ~ Soil * Day
## Value Std.Error DF t-value p-value
## (Intercept) 6.769938 0.03264816 118 207.36050 0e+00
## Soil70 0.307848 0.04617147 22 6.66750 0e+00
## Day 0.015753 0.00268973 118 5.85679 0e+00
## Soil70:Day -0.013615 0.00380385 118 -3.57928 5e-04
## Correlation:
## (Intr) Soil70 Day
## Soil70 -0.707
## Day -0.851 0.602
## Soil70:Day 0.602 -0.851 -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.4222157 -0.2534009 0.1372206 0.6438246 1.5195804
##
## Number of Observations: 144
## Number of Groups: 24
# second approach ---
# full model
lmm0_2 <- lme(log(richness) ~ Soil * Day *Frequency, random=~1|box, data=df)
summary(lmm0_2)
## Linear mixed-effects model fit by REML
## Data: df
## AIC BIC logLik
## -24.5273 26.80925 30.26365
##
## Random effects:
## Formula: ~1 | box
## (Intercept) Residual
## StdDev: 3.785982e-06 0.1420281
##
## Fixed effects: log(richness) ~ Soil * Day * Frequency
## Value Std.Error DF t-value p-value
## (Intercept) 6.670841 0.06380629 112 104.54833 0.0000
## Soil70 0.330921 0.09023572 16 3.66729 0.0021
## Day 0.018025 0.00525670 112 3.42895 0.0008
## Frequencyf2 0.180858 0.09023572 16 2.00428 0.0623
## Frequencyf3 0.134184 0.09023572 16 1.48704 0.1564
## Frequencyf4 0.081348 0.09023572 16 0.90150 0.3807
## Soil70:Day -0.006909 0.00743409 112 -0.92938 0.3547
## Soil70:Frequencyf2 -0.120291 0.12761258 16 -0.94263 0.3599
## Soil70:Frequencyf3 -0.104199 0.12761258 16 -0.81653 0.4262
## Soil70:Frequencyf4 0.132201 0.12761258 16 1.03595 0.3156
## Day:Frequencyf2 -0.005614 0.00743409 112 -0.75514 0.4518
## Day:Frequencyf3 -0.004108 0.00743409 112 -0.55261 0.5816
## Day:Frequencyf4 0.000635 0.00743409 112 0.08538 0.9321
## Soil70:Day:Frequencyf2 -0.005746 0.01051339 112 -0.54658 0.5858
## Soil70:Day:Frequencyf3 -0.002517 0.01051339 112 -0.23944 0.8112
## Soil70:Day:Frequencyf4 -0.018560 0.01051339 112 -1.76536 0.0802
## Correlation:
## (Intr) Soil70 Day Frqnc2 Frqnc3 Frqnc4 Sl70:D S70:F2
## Soil70 -0.707
## Day -0.851 0.602
## Frequencyf2 -0.707 0.500 0.602
## Frequencyf3 -0.707 0.500 0.602 0.500
## Frequencyf4 -0.707 0.500 0.602 0.500 0.500
## Soil70:Day 0.602 -0.851 -0.707 -0.426 -0.426 -0.426
## Soil70:Frequencyf2 0.500 -0.707 -0.426 -0.707 -0.354 -0.354 0.602
## Soil70:Frequencyf3 0.500 -0.707 -0.426 -0.354 -0.707 -0.354 0.602 0.500
## Soil70:Frequencyf4 0.500 -0.707 -0.426 -0.354 -0.354 -0.707 0.602 0.500
## Day:Frequencyf2 0.602 -0.426 -0.707 -0.851 -0.426 -0.426 0.500 0.602
## Day:Frequencyf3 0.602 -0.426 -0.707 -0.426 -0.851 -0.426 0.500 0.301
## Day:Frequencyf4 0.602 -0.426 -0.707 -0.426 -0.426 -0.851 0.500 0.301
## Soil70:Day:Frequencyf2 -0.426 0.602 0.500 0.602 0.301 0.301 -0.707 -0.851
## Soil70:Day:Frequencyf3 -0.426 0.602 0.500 0.301 0.602 0.301 -0.707 -0.426
## Soil70:Day:Frequencyf4 -0.426 0.602 0.500 0.301 0.301 0.602 -0.707 -0.426
## S70:F3 S70:F4 Dy:Fr2 Dy:Fr3 Dy:Fr4 S70:D:F2 S70:D:F3
## Soil70
## Day
## Frequencyf2
## Frequencyf3
## Frequencyf4
## Soil70:Day
## Soil70:Frequencyf2
## Soil70:Frequencyf3
## Soil70:Frequencyf4 0.500
## Day:Frequencyf2 0.301 0.301
## Day:Frequencyf3 0.602 0.301 0.500
## Day:Frequencyf4 0.301 0.602 0.500 0.500
## Soil70:Day:Frequencyf2 -0.426 -0.426 -0.707 -0.354 -0.354
## Soil70:Day:Frequencyf3 -0.851 -0.426 -0.354 -0.707 -0.354 0.500
## Soil70:Day:Frequencyf4 -0.426 -0.851 -0.354 -0.354 -0.707 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.1482469 -0.2683141 0.1338492 0.5505475 1.6809933
##
## Number of Observations: 144
## Number of Groups: 24
# update model
lmm1_2 <- lme(log(richness) ~ Soil * Day, random=~1|box, data=df)
summary(lmm1_2)
## Linear mixed-effects model fit by REML
## Data: df
## AIC BIC logLik
## -106.2061 -88.55627 59.10306
##
## Random effects:
## Formula: ~1 | box
## (Intercept) Residual
## StdDev: 1.134264e-05 0.1453448
##
## Fixed effects: log(richness) ~ Soil * Day
## Value Std.Error DF t-value p-value
## (Intercept) 6.769938 0.03264816 118 207.36050 0e+00
## Soil70 0.307848 0.04617147 22 6.66750 0e+00
## Day 0.015753 0.00268973 118 5.85679 0e+00
## Soil70:Day -0.013615 0.00380385 118 -3.57928 5e-04
## Correlation:
## (Intr) Soil70 Day
## Soil70 -0.707
## Day -0.851 0.602
## Soil70:Day 0.602 -0.851 -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.4222157 -0.2534009 0.1372206 0.6438246 1.5195804
##
## Number of Observations: 144
## Number of Groups: 24
# full model which considers non homogeneouse of variance
gls0 <- gls(log(richness) ~ Soil * Day * Frequency, correlation=corAR1(form= ~Day |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls0)
## Generalized least squares fit by REML
## Model: log(richness) ~ Soil * Day * Frequency
## Data: df
## AIC BIC logLik
## -33.79917 20.38941 35.89958
##
## Correlation Structure: ARMA(1,0)
## Formula: ~Day | box
## Parameter estimate(s):
## Phi1
## 0
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Soil
## Parameter estimates:
## 0 70
## 1.000000 1.530899
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 6.670841 0.04934780 135.18010 0.0000
## Soil70 0.330921 0.09023571 3.66729 0.0004
## Day 0.018025 0.00406553 4.43360 0.0000
## Frequencyf2 0.180858 0.06978833 2.59152 0.0107
## Frequencyf3 0.134184 0.06978833 1.92273 0.0567
## Frequencyf4 0.081348 0.06978833 1.16563 0.2459
## Soil70:Day -0.006909 0.00743409 -0.92938 0.3544
## Soil70:Frequencyf2 -0.120291 0.12761257 -0.94263 0.3476
## Soil70:Frequencyf3 -0.104199 0.12761257 -0.81653 0.4157
## Soil70:Frequencyf4 0.132201 0.12761257 1.03595 0.3022
## Day:Frequencyf2 -0.005614 0.00574953 -0.97639 0.3307
## Day:Frequencyf3 -0.004108 0.00574953 -0.71452 0.4762
## Day:Frequencyf4 0.000635 0.00574953 0.11040 0.9123
## Soil70:Day:Frequencyf2 -0.005746 0.01051339 -0.54658 0.5856
## Soil70:Day:Frequencyf3 -0.002517 0.01051339 -0.23944 0.8111
## Soil70:Day:Frequencyf4 -0.018560 0.01051339 -1.76536 0.0799
##
## Correlation:
## (Intr) Soil70 Day Frqnc2 Frqnc3 Frqnc4 Sl70:D S70:F2
## Soil70 -0.547
## Day -0.851 0.466
## Frequencyf2 -0.707 0.387 0.602
## Frequencyf3 -0.707 0.387 0.602 0.500
## Frequencyf4 -0.707 0.387 0.602 0.500 0.500
## Soil70:Day 0.466 -0.851 -0.547 -0.329 -0.329 -0.329
## Soil70:Frequencyf2 0.387 -0.707 -0.329 -0.547 -0.273 -0.273 0.602
## Soil70:Frequencyf3 0.387 -0.707 -0.329 -0.273 -0.547 -0.273 0.602 0.500
## Soil70:Frequencyf4 0.387 -0.707 -0.329 -0.273 -0.273 -0.547 0.602 0.500
## Day:Frequencyf2 0.602 -0.329 -0.707 -0.851 -0.426 -0.426 0.387 0.466
## Day:Frequencyf3 0.602 -0.329 -0.707 -0.426 -0.851 -0.426 0.387 0.233
## Day:Frequencyf4 0.602 -0.329 -0.707 -0.426 -0.426 -0.851 0.387 0.233
## Soil70:Day:Frequencyf2 -0.329 0.602 0.387 0.466 0.233 0.233 -0.707 -0.851
## Soil70:Day:Frequencyf3 -0.329 0.602 0.387 0.233 0.466 0.233 -0.707 -0.426
## Soil70:Day:Frequencyf4 -0.329 0.602 0.387 0.233 0.233 0.466 -0.707 -0.426
## S70:F3 S70:F4 Dy:Fr2 Dy:Fr3 Dy:Fr4 S70:D:F2 S70:D:F3
## Soil70
## Day
## Frequencyf2
## Frequencyf3
## Frequencyf4
## Soil70:Day
## Soil70:Frequencyf2
## Soil70:Frequencyf3
## Soil70:Frequencyf4 0.500
## Day:Frequencyf2 0.233 0.233
## Day:Frequencyf3 0.466 0.233 0.500
## Day:Frequencyf4 0.233 0.466 0.500 0.500
## Soil70:Day:Frequencyf2 -0.426 -0.426 -0.547 -0.273 -0.273
## Soil70:Day:Frequencyf3 -0.851 -0.426 -0.273 -0.547 -0.273 0.500
## Soil70:Day:Frequencyf4 -0.426 -0.851 -0.273 -0.273 -0.547 0.500 0.500
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.6858047 -0.3073183 0.1130485 0.5706062 1.7771310
##
## Residual standard error: 0.1098446
## Degrees of freedom: 144 total; 128 residual
# full model
glmm1 <- glmer(richness ~ Soil * Day * Frequency + (Day | box), family = poisson, df)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0115456 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
summary(glmm1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ Soil * Day * Frequency + (Day | box)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 3546.2 3602.6 -1754.1 3508.2 125
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -14.2555 -1.2237 0.2777 2.1606 9.9149
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## box (Intercept) 5.994e-03 0.077424
## Day 3.195e-05 0.005652 -0.93
## Number of obs: 144, groups: box, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.681559 0.047209 141.531 < 2e-16 ***
## Soil70 0.345403 0.066308 5.209 1.90e-07 ***
## Day 0.017647 0.003473 5.081 3.75e-07 ***
## Frequencyf2 0.168198 0.066532 2.528 0.01147 *
## Frequencyf3 0.135751 0.066578 2.039 0.04145 *
## Frequencyf4 0.072005 0.066628 1.081 0.27983
## Soil70:Day -0.007814 0.004878 -1.602 0.10923
## Soil70:Frequencyf2 -0.129669 0.093618 -1.385 0.16603
## Soil70:Frequencyf3 -0.110441 0.093651 -1.179 0.23829
## Soil70:Frequencyf4 0.119978 0.093589 1.282 0.19985
## Day:Frequencyf2 -0.004996 0.004896 -1.020 0.30751
## Day:Frequencyf3 -0.004522 0.004899 -0.923 0.35602
## Day:Frequencyf4 0.001281 0.004899 0.262 0.79371
## Soil70:Day:Frequencyf2 -0.003891 0.006893 -0.565 0.57239
## Soil70:Day:Frequencyf3 -0.001878 0.006894 -0.272 0.78532
## Soil70:Day:Frequencyf4 -0.017851 0.006890 -2.591 0.00958 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.0115456 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
# update model
# update model by considering overdispersion
# add a observation-level random effect (sample) if overdispersed Elston et al. (2001).
df$sample <- factor(paste(df$Soil, df$Frequency, df$Day, df$replicates, sep = "_"))
str(df)
## 'data.frame': 144 obs. of 7 variables:
## $ richness : int 828 1075 1101 1100 1108 1124 792 900 727 1074 ...
## $ Soil : Factor w/ 2 levels "0","70": 1 1 1 1 1 1 1 1 1 1 ...
## $ Frequency : Factor w/ 4 levels "f1","f2","f3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : int 12 12 12 16 16 16 2 2 2 20 ...
## $ replicates: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
## $ box : Factor w/ 24 levels "0_f1_A","0_f1_B",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ sample : Factor w/ 144 levels "0_f1_12_A","0_f1_12_B",..: 1 2 3 4 5 6 10 11 12 7 ...
glmm2 <- glmer(richness ~ Soil * Day + (Day | box) + (1 | sample), family = poisson, df)
## boundary (singular) fit: see ?isSingular
# check overdispersion https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-for-overdispersioncomputing-overdispersion-factor
overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
# p < 0.05, overdispersed
overdisp_fun(glmm2)
## chisq ratio rdf p
## 8.46045473 0.06220923 136.00000000 1.00000000
tab_model(glmm2,show.icc=T)
## boundary (singular) fit: see ?isSingular
| richness | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 871.77 | 815.12 – 932.34 | <0.001 |
| Soil [70] | 1.36 | 1.24 – 1.50 | <0.001 |
| Day | 1.02 | 1.01 – 1.02 | <0.001 |
| Soil [70] * Day | 0.99 | 0.98 – 0.99 | <0.001 |
| Random Effects | |||
| σ2 | 0.02 | ||
| τ00 sample | 0.02 | ||
| τ00 box | 0.00 | ||
| τ11 box.Day | 0.00 | ||
| ρ01 box | -1.00 | ||
| N box | 24 | ||
| N sample | 144 | ||
| Observations | 144 | ||
| Marginal R2 / Conditional R2 | 0.389 / NA | ||
summary(glmm2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ Soil * Day + (Day | box) + (1 | sample)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 1883.5 1907.3 -933.8 1867.5 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.24534 -0.05996 0.02397 0.12886 0.28423
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sample (Intercept) 1.837e-02 0.135532
## box (Intercept) 2.333e-03 0.048301
## Day 7.295e-06 0.002701 -1.00
## Number of obs: 144, groups: sample, 144; box, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.770520 0.034277 197.525 < 2e-16 ***
## Soil70 0.308498 0.048363 6.379 1.78e-10 ***
## Day 0.015717 0.002689 5.845 5.07e-09 ***
## Soil70:Day -0.013607 0.003797 -3.584 0.000339 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Soil70 Day
## Soil70 -0.709
## Day -0.863 0.611
## Soil70:Day 0.611 -0.863 -0.708
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(glmm2)
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Soil 1 42.122 42.122 42.122
## Day 1 21.938 21.938 21.938
## Soil:Day 1 12.844 12.844 12.844
glmm的syntax写的是否合适,能否体现时间自相关? summary结果里的Correlation of Fixed Effects很费解啊。。。
# brm0 <- brm(bf(richness ~ Soil + Day + (Day|box), sigma ~ Soil), data = df, autocor = cor_cosy(~Day|box))
# brm1 <- brm(richness ~ Soil * Day, data = df, family = negbinomial("log"))
# summary(brm1)
# update model by removing not significant factors
gls1 <- gls(log(richness) ~ Soil + Day + Soil:Day, correlation=corAR1(form= ~Day |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls1)
## Generalized least squares fit by REML
## Model: log(richness) ~ Soil + Day + Soil:Day
## Data: df
## AIC BIC logLik
## -114.0469 -93.45544 64.02347
##
## Correlation Structure: ARMA(1,0)
## Formula: ~Day | box
## Parameter estimate(s):
## Phi1
## 0
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Soil
## Parameter estimates:
## 0 70
## 1.000000 1.461337
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 6.769938 0.02607474 259.63593 0e+00
## Soil70 0.307848 0.04617147 6.66750 0e+00
## Day 0.015753 0.00214817 7.33328 0e+00
## Soil70:Day -0.013615 0.00380385 -3.57928 5e-04
##
## Correlation:
## (Intr) Soil70 Day
## Soil70 -0.565
## Day -0.851 0.481
## Soil70:Day 0.481 -0.851 -0.565
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.1849151 -0.2226717 0.1280295 0.6379215 1.6976512
##
## Residual standard error: 0.1160809
## Degrees of freedom: 144 total; 140 residual
anova(gls1)
## Denom. DF: 140
## numDF F-value p-value
## (Intercept) 1 382894.3 <.0001
## Soil 1 47.6 <.0001
## Day 1 41.4 <.0001
## Soil:Day 1 12.8 5e-04
plot(gls1)
on day 0, the richness of 0 year soil inundated is \(e^{6.764}\), while the richness of 70 year soil at f1 is \(e^{(6.764+0.331)}\), i.e. \(e^{7.095}\). For the richness of 0 and 70 year soils are:
\(richness_{0yr} = e^{(6.764 + 0.0163*Day)}\); and
\(richness_{70yr} = e^{(7.095 + (0.0163-0.0146)*Day)}\), i.e. \(richness_{70yr} = e^{(7.095 + 0.0018*Day)}\), respectively.
# compare two models
AIC(gls0, gls1)
## Warning in AIC.default(gls0, gls1): models are not all fitted to the same number
## of observations
## df AIC
## gls0 19 -33.79917
## gls1 7 -114.04694
plot(gls1)
# check normality via histogram (1)
hist(gls1$residuals, main="Residual distribution")
## assess the normality of residuals using shapiro wilk test (2)
shapiro.test(gls1$residuals)
##
## Shapiro-Wilk normality test
##
## data: gls1$residuals
## W = 0.84583, p-value = 5.526e-11
plot(gls1$residuals ~ gls1$fitted)
(f <- ggplot(df, aes(x = Day, y = log(richness), shape = Frequency, color = Soil)) +
geom_point(size = 2, alpha=0.9) +
facet_grid(Soil ~ Frequency) +
labs(x = "Days", y = "Richness")+
theme_light())
the Shapiro Wilk test was significant (p < 0.05), indicating residuals are not normal distributed?!
# read dataset
df <- read.csv("alpha_diversity_dispersal.csv", row.names = 1, check.names = 0)
# only keep richness
df <- df[, -c(2:7)]
# remove samples from day 0
df <- subset(df, !grepl('___', rownames(df)))
# select soil from sea water treatment
df <- df[df$Frequency == "f1" | df$Frequency == "f4",]
df$Soil <- factor(df$Soil, levels = c("0", "70"))
df$Water <- factor(df$Water)
df$Frequency <- factor(df$Frequency)
df$box <- factor(paste(df$Soil, df$Frequency, df$replicates, sep = "_"))
str(df)
## 'data.frame': 144 obs. of 7 variables:
## $ richness : int 828 1075 1101 1100 1108 1124 792 900 727 1074 ...
## $ Soil : Factor w/ 2 levels "0","70": 1 1 1 1 1 1 1 1 1 1 ...
## $ Frequency : Factor w/ 2 levels "f1","f4": 1 1 1 1 1 1 1 1 1 1 ...
## $ Water : Factor w/ 2 levels "sea","strl": 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : int 12 12 12 16 16 16 2 2 2 20 ...
## $ replicates: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
## $ box : Factor w/ 12 levels "0_f1_A","0_f1_B",..: 1 2 3 1 2 3 1 2 3 1 ...
head(df)
## richness Soil Frequency Water Day replicates box
## 0_f1_sea_12_A 828 0 f1 sea 12 A 0_f1_A
## 0_f1_sea_12_B 1075 0 f1 sea 12 B 0_f1_B
## 0_f1_sea_12_C 1101 0 f1 sea 12 C 0_f1_C
## 0_f1_sea_16_A 1100 0 f1 sea 16 A 0_f1_A
## 0_f1_sea_16_B 1108 0 f1 sea 16 B 0_f1_B
## 0_f1_sea_16_C 1124 0 f1 sea 16 C 0_f1_C
# histogram
hist(df$richness, main = "histogram of richness")
(f <- ggplot(df, aes(x = Day, y = richness, shape = Frequency, color = Water)) +
geom_point(size = 2, alpha=0.9) +
facet_grid(Soil ~ Frequency) +
labs(x = "Days", y = "Richness")+
theme_light())
(f <- ggplot(df, aes(x = factor(Day), y = richness, shape = Frequency, color = Water)) +
geom_boxplot() +
facet_grid(Soil ~ .) +
geom_point(size = 2, alpha=0.9, position = position_jitterdodge()) +
labs(x = "Days", y = "Richness")+
theme_light())
# full model
gls0 <- gls(log(richness) ~ Soil * Day * Water * Frequency, correlation=corAR1(form= ~1 |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls0)
## Generalized least squares fit by REML
## Model: log(richness) ~ Soil * Day * Water * Frequency
## Data: df
## AIC BIC logLik
## -34.53114 19.65744 36.26557
##
## Correlation Structure: AR(1)
## Formula: ~1 | box
## Parameter estimate(s):
## Phi
## -0.04157485
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Soil
## Parameter estimates:
## 0 70
## 1.0000000 0.7923043
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 6.670819 0.06919295 96.40895 0.0000
## Soil70 0.330203 0.08827854 3.74047 0.0003
## Day 0.018140 0.00579290 3.13139 0.0022
## Waterstrl 0.150218 0.09792631 1.53399 0.1275
## Frequencyf4 0.082946 0.09785360 0.84765 0.3982
## Soil70:Day -0.007070 0.00739076 -0.95662 0.3406
## Soil70:Waterstrl -0.078749 0.12493747 -0.63031 0.5296
## Day:Waterstrl -0.011984 0.00818955 -1.46329 0.1458
## Soil70:Frequencyf4 0.131631 0.12484471 1.05436 0.2937
## Day:Frequencyf4 0.000335 0.00819240 0.04090 0.9674
## Waterstrl:Frequencyf4 -0.215139 0.13848871 -1.55347 0.1228
## Soil70:Day:Waterstrl 0.008981 0.01044849 0.85958 0.3916
## Soil70:Day:Frequencyf4 -0.018248 0.01045212 -1.74583 0.0832
## Soil70:Waterstrl:Frequencyf4 0.032896 0.17668826 0.18618 0.8526
## Day:Waterstrl:Frequencyf4 0.012234 0.01158178 1.05628 0.2928
## Soil70:Day:Waterstrl:Frequencyf4 -0.002956 0.01477640 -0.20003 0.8418
##
## Correlation:
## (Intr) Soil70 Day Wtrstr Frqnc4 Sl70:D
## Soil70 -0.784
## Day -0.865 0.678
## Waterstrl -0.708 0.555 0.612
## Frequencyf4 -0.707 0.554 0.612 0.500
## Soil70:Day 0.678 -0.865 -0.784 -0.480 -0.480
## Soil70:Waterstrl 0.555 -0.708 -0.480 -0.784 -0.392 0.612
## Day:Waterstrl 0.611 -0.479 -0.707 -0.865 -0.432 0.554
## Soil70:Frequencyf4 0.554 -0.707 -0.480 -0.392 -0.784 0.612
## Day:Frequencyf4 0.612 -0.480 -0.707 -0.433 -0.865 0.554
## Waterstrl:Frequencyf4 0.500 -0.392 -0.433 -0.707 -0.708 0.339
## Soil70:Day:Waterstrl -0.479 0.611 0.554 0.678 0.339 -0.707
## Soil70:Day:Frequencyf4 -0.480 0.612 0.554 0.339 0.678 -0.707
## Soil70:Waterstrl:Frequencyf4 -0.392 0.500 0.339 0.554 0.555 -0.433
## Day:Waterstrl:Frequencyf4 -0.432 0.339 0.500 0.611 0.611 -0.392
## Soil70:Day:Waterstrl:Frequencyf4 0.339 -0.432 -0.392 -0.479 -0.479 0.500
## Sl70:W Dy:Wtr S70:F4 Dy:Fr4 Wtr:F4 Sl70:D:W
## Soil70
## Day
## Waterstrl
## Frequencyf4
## Soil70:Day
## Soil70:Waterstrl
## Day:Waterstrl 0.678
## Soil70:Frequencyf4 0.500 0.339
## Day:Frequencyf4 0.339 0.500 0.678
## Waterstrl:Frequencyf4 0.554 0.611 0.555 0.612
## Soil70:Day:Waterstrl -0.865 -0.784 -0.432 -0.392 -0.479
## Soil70:Day:Frequencyf4 -0.433 -0.392 -0.865 -0.784 -0.480 0.500
## Soil70:Waterstrl:Frequencyf4 -0.707 -0.479 -0.708 -0.480 -0.784 0.611
## Day:Waterstrl:Frequencyf4 -0.479 -0.707 -0.479 -0.707 -0.865 0.554
## Soil70:Day:Waterstrl:Frequencyf4 0.611 0.554 0.611 0.554 0.678 -0.707
## S70:D:F S70:W: D:W:F4
## Soil70
## Day
## Waterstrl
## Frequencyf4
## Soil70:Day
## Soil70:Waterstrl
## Day:Waterstrl
## Soil70:Frequencyf4
## Day:Frequencyf4
## Waterstrl:Frequencyf4
## Soil70:Day:Waterstrl
## Soil70:Day:Frequencyf4
## Soil70:Waterstrl:Frequencyf4 0.612
## Day:Waterstrl:Frequencyf4 0.554 0.678
## Soil70:Day:Waterstrl:Frequencyf4 -0.707 -0.865 -0.784
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.6712668 -0.4130164 0.1313822 0.6485909 1.4405870
##
## Residual standard error: 0.152314
## Degrees of freedom: 144 total; 128 residual
# update model by removing not significant factors
gls1 <- gls(log(richness) ~ Soil + Day + Soil:Day, correlation=corAR1(form= ~1 |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls1)
## Generalized least squares fit by REML
## Model: log(richness) ~ Soil + Day + Soil:Day
## Data: df
## AIC BIC logLik
## -117.3097 -96.71821 65.65486
##
## Correlation Structure: AR(1)
## Formula: ~1 | box
## Parameter estimate(s):
## Phi
## -0.07361147
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Soil
## Parameter estimates:
## 0 70
## 1.0000000 0.8203663
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 6.733888 0.03512164 191.73046 0.0000
## Soil70 0.364713 0.04542792 8.02838 0.0000
## Day 0.015352 0.00297846 5.15439 0.0000
## Soil70:Day -0.012466 0.00385247 -3.23582 0.0015
##
## Correlation:
## (Intr) Soil70 Day
## Soil70 -0.773
## Day -0.877 0.678
## Soil70:Day 0.678 -0.877 -0.773
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.9921583 -0.4655684 0.1743520 0.6703851 1.6376582
##
## Residual standard error: 0.1534817
## Degrees of freedom: 144 total; 140 residual
anova(gls1)
## Denom. DF: 140
## numDF F-value p-value
## (Intercept) 1 430638.5 <.0001
## Soil 1 116.4 <.0001
## Day 1 17.5 0.0001
## Soil:Day 1 10.5 0.0015
on day 0, the richness of 0 year soil inundated is \(e^{6.733}\), while the richness of 70 year soil at f1 is \(e^{(6.733+0.365)}\), i.e. \(e^{7.098}\). For the richness of 0 and 70 year soils are:
\(richness_{0yr} = e^{(6.733 + 0.0154*Day)}\); and
\(richness_{70yr} = e^{(7.098 + (0.0154-0.0124)*Day)}\), i.e. \(richness_{70yr} = e^{(7.098 + 0.0030*Day)}\), respectively.
# compare two models
AIC(gls0, gls1)
## Warning in AIC.default(gls0, gls1): models are not all fitted to the same number
## of observations
## df AIC
## gls0 19 -34.53114
## gls1 7 -117.30971
plot(gls1)
# check normality via histogram (1)
hist(gls1$residuals, main="Residual distribution")
## assess the normality of residuals using shapiro wilk test (2)
shapiro.test(gls1$residuals)
##
## Shapiro-Wilk normality test
##
## data: gls1$residuals
## W = 0.90314, p-value = 3.324e-08
plot(gls1$residuals ~ gls1$fitted)
(f <- ggplot(df, aes(x = Day, y = log(richness), shape = Frequency, color = Soil)) +
geom_point(size = 2, alpha=0.9) +
facet_grid(Soil ~ Frequency) +
labs(x = "Days", y = "Richness")+
theme_light())
the Shapiro Wilk test was significant (p < 0.05), indicating residuals are not normal distributed?!